Real-Time Diagrammatic Monte Carlo for Nonequilibrium Quantum Transport 
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We propose a novel approach to nonequilibrium real-time dynamics of quantum impurities models 
coupled to biased non-interacting leads, such as those relevant to quantum transport in nanoscale 
molecular devices. The method is based on a Diagrammatic Monte Carlo sampling of the real-time 
perturbation theory along the Keldysh contour. We benchmark the method on a non-interacting res- 
onant level model and, as a first non-trivial application, we study zero temperature non-equilibrium 
transport through a vibrating molecule. 
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Introduction. Recent advances in nanotechnology have 
made it possible to contact microscopic quantum objects, 
like artificial atoms (quantum dots), molecules or quan- 
tum wires, with metallic electrodes, opening the route 
towards very promising nanoelectronical devices [l], Q. 
Following the first discovery of the Kondo effect in quan- 
tum dots [3[ , a lot of challenging experiments have been 
carried out contacting^ e.g., single-molecules or cou- 
pled quantum dots [6| to metallic leads. In the mean- 
time, these impressive experimental developments have 
raised a lot of novel and interesting physical questions, 
so that electronic transport through nano-systems has 
become one of the frontier fields in condensed matter 
physics. Indeed, the possibility contact microscopic ob- 
jects is particulary intriguing as it allows to study quan- 
tum transport in a regime where the tunneling rate be- 
comes comparable or even smaller than other energy 
scales, like the electron-electron repulsion or the energy 
of atomic displacements .which may lead to interesting 
non-equilibrium effects. 0, @] 

These experimental progresses urgently ask for de- 
veloping efficient non-perturbative theoretical tools to 
treat out-of-equilibrium phenomena. The simplest way 
to model nonequilibrium transport in nanodevices is 
through a quantum impurity model, namely a set of dis- 
crete levels a (with creation operators c' a(X , a being the 
spin), mimicking a quantum dot or a molecule, coupled to 
two baths of non interacting electrons (with creation op- 
erators fl aa ), labeled by some quantum number k, which 
account for the metallic leads. The leads (a = L,R) are 
kept at different chemical potentials /jl — /zr = eV. As 
a consequence the general Hamiltonian may read 



^ — ^ £ ka ftcrg fkcra + ^loc \c\g ; c aa] 

a—L.R k a 

E (Vka a fL a C at ,-rh.c). (1) 



The local Hamiltonian 7ii oc accounts for all the physics 
on the quantum impurity, including possibly vibrational 



degrees of freedoms, and, because of the discrete set of 
levels, could in principle be diagonalized exactly. How- 
ever, the hybridization to the reservoirs makes the prob- 
lem untreatable unless in simple cases. Furthermore, the 
finite bias, which drives the system out of equilibrium, 
rules out the possibility to apply all the powerful tools 
developed during last decades, like Numerical Renormal- 
ization Group (NRG) 0] and the recent Diagrammatic 
Monte Carlo method (DiagMC) [13, El- 
Standard approaches to nonequilibrium are usually based 
on Keldysh perturbation theory [12j, which is analyti- 
cally feasible only to lowest orders. In order to access to 
the most interesting intermediate-coupling regime several 
numerical methods has been proposed 1J, 
In this Letter we present a novel real-time DiagMC ap- 
proach to nonequilibrium transport in quantum impurity 
models. The method is based on a stochastic sampling 
of the Keldysh diagrams generated by the perturbative 
expansion in the tunnelling. It does not require any dis- 
cretization of the time evolution, hence it provides a very 
accurate description of the dynamics. We benchmark the 
method on a non-interacting resonant level and then, as 
a first non-trivial application, we compute the inelastic 
tunneling spectrum in the resonant level coupled to a lo- 
cal vibrational mode. 

Formulation. To set up the method, we consider an ini- 
tially decoupled system, made by the isolated impurity 
and the two leads, each assumed to be at equilibrium 
with its own reservoir at chemical potential fi a . At time 
t > 0, we switch-on the hybridization in |T]) and let the 
system evolve with the full hamiltonian TL. Given an 
initial density matrix, po, we want to compute average 
values of physical operators evolved in real-time from 
to t. The final goal is to succeed following the dynamics 
till the steady state, so to compute experimentally rele- 
vant quantities like the differential conductance dl /dV. 

Real-time quantum dynamics can be generally rep- 
resented as evolution along the so-called Keldysh con- 
tour /C, plotted in Fig. [TJ made of two branches winding 
around the real-time axis from to t (lower branch) and 
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FIG. 1: a) Keldysh contour /C used in real-time diagMC. 
Under time-ordering Tk. creation/annihilation operators are 
displaced along the contour time-ordered as shown by the ar- 
rows, b) An example of a configuration with k = 3 segments, 
each starting at tf and ending at t\ , i — 1,2, 3. Here blue/red 
dots stands for annihilation/creation operators of the initially 
occupied level, while red segments indicate how the vertices 
are connected by the hybridization functions Aje(t e , t 3 ) in the 
particular configuration shown. 



back from t to (upper branch). In this perspective the 
average value of any operator O can be written as 

(O (t)) = Tr [ Po T K (e^h J , ( 2 ) 

where the trace is over the lead and impurity degrees of 



freedom and Tjq denotes time-ordering along the Keldysh 
contour. The main idea of the approach is to expand 
the evolution operator @ in powers of the hybridization 
and trace-out exactly the leads degrees of freedom. The 
expansion obtained in that way looks as the natural gen- 
eralization of the diagrammatic expansion of Ref. 11| to 
the /C contour, which is required to deal with nonequi- 
librium effects. The resulting diagrams are then sam- 
pled with a Monte Carlo algorithm that we shall dis- 
cuss in detail. To show how the method works, we start 
by considering a spinless biased resonant level model 
(RLM), namely a single fermionic energy level Ed driven 
out of equilibrium by an applied bias eV — [II — pn 
between the two leads. Given an initial density matrix 
Po = Pleads ® Pimp, with p leads describing the two un- 
coupled leads each at equilibrium with its own chemical 
potential and pi mp = c'c — the impurity is initially occu- 
pied - we are interested in computing, for example, the 
real-time evolution of the RLM population n = c'c. To 
this extent, we expand the Keldysh evolution operator in 
powers of the hybrization between the bath and the im- 
purity and trace out the lead degrees of freedom, which 
are non-interacting. The resulting expansion reads 
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Here T>u is the outcome of tracing the lead degrees of 
freedom and can be expressed in closed form as the de- 
terminant of a k x k matrix Mr 1 , 

V(tt,...,tt\t{,...,t%)=det (M- 1 ), (4) 

whose entries are the Keldysh hybridization functions 

M^=i^ K {tt,tf) a(f?,tj), (5) 

which we define as 

A K = £ \V ka f(T K (f ka (O/L C s ))> • (6) 

ka 

Here we adopt the standard definition of the Keldysh 
Green's functions [16(, namely we consider t s , t e as living 
on the contour IC. We note the additional sign s (t e ,t s ), 
which is negative when the two times are on opposite 
branches and positive otherwise. While the determinant 
T>h properly accounts for the effects associated to the 



leads, the function Ck involves only the impurity degrees 
of freedom and can be generally written as a trace over 
the initial impurity density matrix, namely 

C h = TrL imp T K (c\tt)c(t s k ) . . . c\t\)c{t\)n{t)) j. 

(7) 

The expansion thus obtained admits a natural represen- 
tation in terms of a collection of k segments, t € [if, tf\ , or 
equivalently k — 1 anti-segments, t 6 [tf,tf +1 ], properly 
ordered along the contour IC and connected in all possible 
ways by the Keldysh hybridization functions A/c(t e ,t s ). 
An example of such a configuration for k = 3 is plotted 
in fig. [T] It is worth noticing that a similar expansion 
can be carried out also for other quantities like e.g. the 
average current (I a (t)) or the noise S(t) = (I a (t)I a (0)), 
resulting into a very general approach. 
Algorithm. As usual in DiagMC [l(| E3 we view the 
perturbative expansion as a sum over configurations C, 
i.e. diagrams with k segments placed along the contour 
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K,. In order to get an efficient sampling scheme three 
basic updates are implemented: adding/removing a seg- 
ment, adding/removing an antisegment or shifting a seg- 
ment end-point. We accept/reject a new configuration 
according to detailed balance prescription. In the actual 
simulation, we store and update the matrix M. defined 
in ([5]), which is the only quantity required to compute 
Metropolis acceptance ratios [10 ]. 

Benchmark. We benchmark the method in the bi- 
ased spinless RLM, which can be exactly solved hence 
being the simplest test for a real-time diagMC calcula- 
tion. In fact, the main issue when sampling real-time 
quantum dynamics is the so called sign-problem, namely 
the exponential increase of relative errors in the Monte 
Carlo estimate of any observable in the infinite-size (or 
zero temperature) limit. This is well known in thermal- 
equilibrium (imaginary-time) Monte Carlo simulations 
of fermionic systems. In this perspective, diagrammatic 
Monte Carlo has been proved to be a generic sign-tolerant 
approach which can deal with sign-alterning series and 
even take advantage from them [171 ] . How far this ap- 
proach can be pushed, in particular in studying real-time 
dynamics, is the main scope of this work. To this extent, 
we compute the occupation number (n(t)) as a function 
of time t for different values of the level position Ed both 
in equilibrium, eV = 0, as well as out of equilibrium, 
eV 7^ 0, and compare the exact results with diagMC 
data. In this simple case, we expect that a single energy 
scale, namely the level broadening T = 7r^ fe |Vfc| 2 (5(£fe), 
controls the approach to the steady-state. As can be seen 
from Fig. [2j this is indeed confirmed by diagMC calcula- 
tion which perfectly matches the exact solution. Indeed, 
we are able to resolve both the short-time decay from the 
initial configuration as well as the approach to the steady 
state. We note that a finite applied bias, eV / 0, cuts 
off the Keldysh evolution operator (JSJ) (the steady state 
is reached earlier than at eV = 0), as pointed out by 



Ref. |18|, thus making the expansion more convergent. 



As mentioned previously, within the present approach 
one can easily measure the current flowing thorugh the 
impurity I(t) = \{Ih{t) — I nit)) on a fine real-time grid 
in a very efficient way, which is also shown in Fig. [5] It 
is worth mentioning that within our Keldysh diagMC we 
are able to reach a true non- equilibrium steady state with 
a finite value of the current, due to the infinite size of the 
bath. Dissipation occurs entirely within the fermionic 
reservoir and we do not need to include any ficticious 
bosonic bath to reach a steady state. 
Nonequilibrium Transport through a single-molecule. As 
a first non trivial application we consider a simple model 
of a molecular conductor, namely a spinless fermionic 
level coupled to Holstein phonon. The local Hamiltonian 
reads 
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FIG. 2: Zero-temperature real-time dynamics for different 
bias values eV of the dot population, {n(t)), from an ini- 
tially occupied dot (left panel) and of the current, (I(t)} (right 
panel) of the resonant level model. DiagMC results (dots) are 
compared with the exact solution (full line). We take ed = T 
and consider a flat density of states in the leads with an half- 
bandwidth tor. 



where ujq is the phonon frequency (with displacement x 
and its conjugate variable p), Ed is the energy of the level 
and n its occupancy and finally g the electron-phonon 
coupling. Our Keldysh diagMC can be naturally ex- 
tendend to include local phonons, the only difference 
appears in the trace over local degrees of freedom 0, 
which now involves fermionic operators evolved accord- 
ing to the hamiltonian ([8]) for the electron-phonon sub- 
system. This trace can be evaluated analytically by ob- 
serving that the local Hamiltonians with different level 
occupancy n — 0, 1 are related one to the other by a uni- 
tary transformation, 7i; oc (0) + ed = J7 T W/ c(l)C/, with 
U — exp (igp/ujo). It follows that the bosonic contribu- 
tion to the local trace reduces to the following bosonic 
correlation-function 



C{ h = Tr[p ph Ui(t%)U(t° k )tf(t%_ 1 ) 



,U(t{) 



(9) 



Hi oc (n) 



w 



(x 2 +p 2 )+gx(n- -)+ £d (n- -) (8) 



which can be easily evaluated analytically for most com- 
mon initial phonon density matrices p p h, which we as- 
sume the equlibrium distribution at zero temperature. 
In ([9|) U(t) and U'(t) are the unitary operators evolved 
with 7i/ oc (l), and we have assumed the level initially oc- 
cupied. Therefore local vibrational degrees of freedom 
can be included in our sampling scheme at any order in 
the coupling constant g, allowing to extract by Keldysh 
DiagMC non-perturbative results in the electron-phonon 
coupling. 

The coupling to molecular vibrations is known to sig- 
nificantly affect inelastic electron tunneling. [3] When 
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FIG. 3: Zero-temperature differential conductance dl /dV (in 
unit of e 2 /h) for bias values eV around wo = 2.0r. Electron- 
phonon coupling is g — 0.5tJo- We consider two different 
values of the fermionic level position Ed in order to reproduce 
the crossover from reducted (step-down) to enhanced (step- 
up) conductance. 



the bias hits a vibrational frequency, the differential con- 
ductance dl/dV changes sharply. Experimentally it is 
observed that dl/dV increases in the tunneling regime, 
but decreases in the opposite case of a high transmis- 
sion barrier. Although simple physical arguments can be 
invoked EH, [13 to expla in this phenomenon, theoret- 
ical calculations within the Keldysh formalism have so 
bar been restricted to the lowest orders in the eiectron- 
HS, 



phonon coupling. [23j,|24, 25] In the simple resonant level 
model that we are considering, the perturbative calcula- 
tions predict that dl/dV, at bias eV = ujo, should de- 
crease if the zero-bias conductance G > 0.5, in units of 
the unitary value (in our spinless case e 2 /h), otherwise 



should increase. [23|, |24j, |25| Evidences in favor of this re- 
sults have been very recently found by Tal and coworkers 
with H 2 molecules bridging a Pt break-junction. Since 
the Keldysh diagMC is, as mentioned, non perturbative 
in the electron-phonon coupling, it offers the possibility 
to verify the above theoretical predictions. We model the 
two regimes of G ^ 0.5 by two different values of the level 
position Ed = 1 and 3, the former closer to resonance than 
the latter, and compute directly the differential conduc- 
tance dl/dV. As can be seen from Fig. |3]a step-down or 
step-up features do appear around the threshold for vi- 
bronic excitations, eV ~ u>q, when the zero bias conduc- 
tance is greater or lower than 0.5, respectively, in agree- 
ment with perturbative results. U 3, [25| We also note 
that the step is not as abrupt as found in perturbation 
theory [23, 24], 25 1, likely signaling a significant phonon 
damping. 

In summary, we have introduced a novel non pertur- 



bative approach to nonequilibrium quantum transport in 
nanoscopic conductors. The method is based on a Di- 
agrammatic Monte Carlo sampling of the real-time per- 
turbation theory in the impurity-leads hybridization, per- 
formed along the Keldysh contour required to treat non 
equilibrium effects. In spite of the oscillating nature of 
the real-time quantum evolution, we are able to follow 
the dynamics starting from an arbitrary initial prepara- 
tion up to the steady state. This is primarily due to 
the combined effects of infinite leads and of the applied 
bias, which cut-off the Keldysh evolution operator, but 
also to the capability of the algorithm to cope with sign- 
problem. As a first application we have studied zero tem- 
perature non linear transport through a simple model 
of a molecular conductor. Being a completely general 
method, it can be in principle used to study any discrete 
quantum system bridging two non-interacting conducting 
leads, providing a new tool to study quantum transport 
in nanoscopic devices. 
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